%pylab inline
import matplotlib.pyplot as plt
from ipywidgets import *
#from scipy import special
#import scipy.special as sp
from scipy.special import jn,fresnel
%pylab is deprecated, use %matplotlib inline and import the required libraries. Populating the interactive namespace from numpy and matplotlib
Talbot effect, near field Fresnel diffraction¶
Talbot carpet, numerical model: ¶
$U(x,z)= \sum_{n=-N}^{N}\sum_{m=-N_s}^{N_s} \frac{e^{i k r_{nm}}}{\sqrt{r_{nm}}}$,¶
where $r_{nm} = \sqrt{z^2+{\left(x-n a+m \frac{d}{2 N_s}\right)}^2}$ . Here $2N+1$ is the number of periods and $2N_s+1$ is the number of test points in one slit.
Theory ¶
The Helmholtz equation: $(\Delta +k^2)U(x,z) =0 $.¶
grating: $U(x,z=0)= T(x) = \sum_{n=-\infty}^{\infty} c_n \, e^{i 2\pi n \frac{x}{a}} \equiv \sum_{n=-\infty}^{\infty} c_n \, e^{i q_n x}$, where $q_n = \frac{2\pi n}{a}$.¶
The total solution satisfying the Helmholtz equation:¶
$U(x,z)= \sum_{n=-\infty}^{\infty} c_n \, e^{i q_n x}e^{i\sqrt{k^2-q_n^2}z}$¶
In paraxial approximation $q_n \ll k$, thus $\sqrt{k^2-q_n^2} \approx k-\frac{q_n^2}{2k}$ and we have¶
$U(x,z) = \sum_{n=-\infty}^{\infty} c_n , e^{i q_n x}e^{i\sqrt{k^2-q_n^2}z} \approx e^{i k z} \sum_{n=-\infty}^{\infty} c_n , e^{i q_n x} e^{-i\frac{q_n^2}{2k} z} \¶
= e^{i k z} \sum_{n=-\infty}^{\infty} c_n , e^{i 2\pi n \frac{x}{a}}e^{-i \pi n^2 \frac{z}{z_T}},$
where $z_T = \frac{a^2}{\lambda}$ is the Talbot period.¶
Note that we can omit the factor $e^{i k z}$ since it drops out in the intensity.
If $z = 2m z_T$ ($m$ is an integer) then $U(x,z+2m z_T) = \sum_{n=-\infty}^{\infty} c_n ,¶
e^{i 2\pi n \frac{x}{a}}e^{-i \pi n^2 \frac{z}{z_T}}e^{-i 2\pi n^2 m} = \sum_{n=-\infty}^{\infty} c_n , e^{i 2\pi n \frac{x}{a}}e^{-i \pi n^2 \frac{z}{z_T}}{\left(e^{-i 2\pi m}\right)}^{n^2} = U(x,z) $.
Note that we omited the factor $e^{i k z}$ since it drops out in the intensity.
This is the revival of the intensity pattern at $z = 2m z_T$.¶
Moreover, it can be shown:¶
\begin{align} U(x,z+z_T) &= U(x+\frac{a}{2},z), \\ {|U(x,z+z_T)|}^2 &= {|U(x+\frac{a}{2},z)|}^2 , \\ {|U(x,z)|}^2 &= {|U(x+\frac{a}{2},z+z_T)|}^2,\\ {|U(x,z)|}^2 &= {|U(x,z+2 z_T)|}^2 . \end{align}
The intensity is periodic with the Talbot period, apart from a shift of $a/2$ between successive images.
Papers:¶
M V Berry and J Goldberg: Renormalisation of curlicues?, Nonlinearity 1, 001 -26 (1988)¶
Berry, M V, Marzoli, I and Schleich, W, 2001 ‘Quantum carpets, carpets of light’ Physics World (June), 39-44 (2001)¶
Berry, M V and Bodenschatz, E, 1999 ‘Caustics, multiply-reconstructed by Talbot interference’, J.Mod.Opt 46 349 – 365¶
def psi2(x,z,param):
# lam ---> hullamhosz, k=2*pi/lam ---> hullamszam,
# a ---> racsallando,
# 2*Nmax ---> racsok szama
# d --> width of one slit
# Ns --> 2*Ns+1 number of test points in one slit
lam, a, d, Nmax, Ns = param
k = 2*pi/lam
sum = 0
for n in range(-Nmax,Nmax+1):
for m in range(-Ns,Ns+1):
rnm = sqrt(z**2+(x-n*a+m*d/2/Ns)**2)
sum +=exp(1j*k*rnm)/sqrt(rnm)
psi = sum
return abs(psi)**2
def spiral(x,z,param):
# lam ---> hullamhosz, k=2*pi/lam ---> hullamszam,
# a ---> racsallando,
# 2*Nmax ---> racsok szama
# d --> width of one slit
# Ns --> 2*Ns+1 number of test points in one slit
lam, a, d, Nmax, Ns = param
k = 2*pi/lam
vec = []
for n in range(-Nmax,Nmax+1):
for m in range(-Ns,Ns+1):
rnm = sqrt(z**2+(x-n*a+m*d/2/Ns)**2)
vec.append(exp(1j*k*rnm)/sqrt(rnm))
return vec
lam,a, d, Nmax, Nw =(0.5,0.1,0.01, 20,1)
zT=a**2/lam ## Talbot-periodus z-ben
print('Talbot period, z_T/a = ' ,zT/a, "\n")
param= [lam,a,d, Nmax, Nw]
psi2(0.2,2,param)
Talbot period, z_T/a = 0.20000000000000004
484.38866410569324
.../okt/Optika_SpR/Talbot/Talbot_with_water_waves_exp_theor.pdf¶
############ near field Fresnel diffraction ############
lam, a, d, Nmax, Nw =(0.01,1,0.1,40,10)
#lam, a, d, Nmax, Nw =(0.01,1,0.35,40,10)
zT=1*a**2/lam ## Talbot-periodus z-ben
xmax,zmax = (1.5*a,6.2*zT)
Np=150
################################################
print('z_T, Talbot period along z = ' ,zT, "\n")
print('Talbot period, z_T/a = ' ,zT/a, "\n")
param= [lam,a,d, Nmax, Nw]
#mintavételezési pontok legyártása
xx,zz = meshgrid(linspace(-xmax,xmax,2*Np),linspace(zT*1.25,zmax,2*Np))
fv= psi2(xx,zz,param)
fvmax=max(flatten(fv))
print('|psi_max|^2 =',fvmax)
### plotting
fig, ax = plt.subplots(1,figsize=(8,6))
res=contourf(xx/a,zz/zT,fv/fvmax,levels=linspace(0.,1,20),cmap='jet') # cmap='rainbow', cmap='viridis', cmap='Reds' cmap='gray'
colorbar(res,fraction=0.049, pad=0.07)
ax.set_xlabel(r'$x/a$',fontsize=12)
ax.set_ylabel(r'$z/z_T$',fontsize=12);
ax.set_title('Near field pattern: Fresnel diffraction')
#ax.set_aspect('equal')
z_T, Talbot period along z = 100.0 Talbot period, z_T/a = 100.0 |psi_max|^2 = 548.8409740515345
Text(0.5, 1.0, 'Near field pattern: Fresnel diffraction')
############ near field Fresnel diffraction ############
#lam, a, d, Nmax, Nw =(0.01,1,0.1,40,10)
lam, a, d, Nmax, Nw =(0.01,1,0.35,40,10)
zT=1*a**2/lam ## Talbot-periodus z-ben
xmax,zmax = (1.5*a,6.2*zT)
Np=150
################################################
print('z_T, Talbot period along z = ' ,zT, "\n")
print('Talbot period, z_T/a = ' ,zT/a, "\n")
param= [lam,a,d, Nmax, Nw]
#mintavételezési pontok legyártása
xx,zz = meshgrid(linspace(-xmax,xmax,2*Np),linspace(zT*1.5,zmax,2*Np))
fv= psi2(xx,zz,param)
fvmax=max(flatten(fv))
print('|psi_max|^2 =',fvmax)
### plotting
fig, ax = plt.subplots(1,figsize=(8,6))
res=contourf(xx/a,zz/zT,fv/fvmax,levels=linspace(0.,1,20),cmap='jet') # cmap='rainbow', cmap='viridis', cmap='Reds' cmap='gray'
colorbar(res,fraction=0.049, pad=0.07)
ax.set_xlabel(r'$x/a$',fontsize=12)
ax.set_ylabel(r'$z/z_T$',fontsize=12);
ax.set_title('Near field pattern: Fresnel diffraction')
#ax.set_aspect('equal')
z_T, Talbot period along z = 100.0 Talbot period, z_T/a = 100.0 |psi_max|^2 = 65.20940140379872
Text(0.5, 1.0, 'Near field pattern: Fresnel diffraction')
############ Fraunhofer diffraction ############
lam, a, d, Nmax, Nw =(0.01,1,0.25,40,10)
zT=1*a**2/lam ## Talbot-periodus z-ben
xmax,zmax = (1.5*a,27.*zT)
Np=75
################################################
print('z_T, Talbot period along z = ' ,zT, "\n")
print('Talbot period, z_T/a = ' ,zT/a, "\n")
param= [lam,a,d, Nmax, Nw]
#mintavételezési pontok legyártása
xx,zz = meshgrid(linspace(-xmax,xmax,2*Np),linspace(zT*1.5,zmax,2*Np))
fv= psi2(xx,zz,param)
fvmax=max(flatten(fv))
print('|psi_max|^2 =',fvmax)
### plotting
fig, ax = plt.subplots(1,figsize=(8,6))
res=contourf(xx/a,zz/zT,fv/fvmax,levels=linspace(0.,1,20),cmap='jet') # cmap='rainbow', cmap='viridis', cmap='Reds' cmap='gray'
colorbar(res,fraction=0.049, pad=0.07)
ax.set_title('Far field pattern: Fraunhofer diffraction')
ax.set_xlabel(r'$x/a$',fontsize=12)
ax.set_ylabel(r'$z/z_T$',fontsize=12);
#ax.set_aspect('equal')
z_T, Talbot period along z = 100.0 Talbot period, z_T/a = 100.0 |psi_max|^2 = 97.71345703255959
lam,a, d, Nmax, Ns =(0.025,1,0.25, 10, 5)
zT=a**2/lam ## Talbot-periodus z-ben
Xp,Zp=(0.0*a,2.0*zT)
#lam,a, d, Nmax, Ns =(0.025,1,0.25, 10, 5)
#zT=a**2/lam ## Talbot-periodus z-ben
#Xp,Zp=(0.0*a,2.0*zT)
# -------------------
#lam,a, d, Nmax, Ns =(0.025,1,0.25, 10, 20)
#zT=a**2/lam ## Talbot-periodus z-ben
#Xp,Zp=(0.*a,2.0*zT)
# -------------------
#lam,a, d, Nmax, Ns =(0.025,1,0.25, 10, 20)
#zT=a**2/lam ## Talbot-periodus z-ben
#Xp,Zp=(0.2*a,2.0*zT)
# -------------------
param= [lam,a,d, Nmax, Ns]
Ap= spiral(Xp,Zp,param)
#Ap, real(Ap), imag(Ap)
lanc=cumsum(Ap,axis=0)
vx,vy=(real(lanc),imag(lanc))
eredo = array([ vx[-1]-vx[0],vy[-1]-vy[0]])
fig, ax = plt.subplots(1, 1,figsize=(8,8))
ax.quiver(vx[:-1], vy[:-1], vx[1:]-vx[:-1], vy[1:]-vy[:-1], color='r',
scale_units='xy', angles='xy', scale=1, width=0.0045)
# ez is kell helyes plottolashoz:
#plot(x[0],y[0], szin, lw=0.01)
ax.quiver(vx[0],vy[0],vx[-1]-vx[0],vy[-1],vy[0], scale_units='xy', angles='xy', scale=1, width=0.005,color='k')
ax.set_aspect('equal')
fact=1.1
ax.set_xlim(fact*min(vx),fact*max(vx))
ax.set_ylim(fact*min(vy),fact*max(vy))
ax.axis('off')
print('Xp= ',Xp/a,'*a \nZp=',Zp/zT,'* zT')
print('Talbot period, z_T/a = ' ,zT, "\n")
print('eredo = ', eredo, norm(eredo))
Xp= 0.0 *a Zp= 2.0 * zT Talbot period, z_T/a = 40.0 eredo = [4.5721166 5.05812758] 6.818277260496455
def rajz_spiral(Xp,Zpp,lam,a,d,Nmax,Ns):
#lam,a, d, Nmax, Ns =(0.025,1,0.25, 10, 5)
# Xp in units of a
# Zpp in units of zT
zT=a**2/lam ## Talbot-periodus z-ben
param= [lam,a,d, Nmax, Ns]
Ap= spiral(Xp*a,Zpp*zT,param)
lanc=cumsum(Ap,axis=0)
vx,vy=(real(lanc),imag(lanc))
eredo = array([ vx[-1]-vx[0],vy[-1]-vy[0]])
fig, ax = plt.subplots(1, 1,figsize=(8,8))
ax.quiver(vx[:-1], vy[:-1], vx[1:]-vx[:-1], vy[1:]-vy[:-1], color='r',
scale_units='xy', angles='xy', scale=1, width=0.0045)
# ez is kell helyes plottolashoz:
#plot(x[0],y[0], szin, lw=0.01)
ax.quiver(vx[0],vy[0],vx[-1]-vx[0],vy[-1],vy[0], scale_units='xy', angles='xy', scale=1, width=0.005,color='k')
ax.set_aspect('equal')
fact=1.1
ax.set_xlim(fact*min(vx),fact*max(vx))
ax.set_ylim(fact*min(vy),fact*max(vy))
ax.axis('off')
print('Xp= ',Xp/a,'*a \nZp=',Zpp,'* zT')
print('Talbot period, z_T/a = ' ,zT, "\n")
print('eredo = ', eredo, norm(eredo))
return
lam,a, d, Nmax, Ns =(0.025,1,0.25, 10, 5)
#zT=a**2/lam ## Talbot-periodus z-ben
Xp,Zpp=(0.0*a,2.0)
rajz_spiral(Xp,Zpp,lam,a,d,Nmax,Ns);
Xp= 0.0 *a Zp= 2.0 * zT Talbot period, z_T/a = 40.0 eredo = [4.5721166 5.05812758] 6.818277260496455
interact(rajz_spiral,
Xp = FloatSlider(min=0,max=2, step=0.1,value=0.0,description=r'$\frac{X_p}{a}$'),
Zpp = FloatSlider(min=1,max=7 ,step=0.5,value=2.0,description=r'$Z_p/z_T$'),
Nmax = IntSlider(min=10,max=60 ,step=5,value=10,description='Nmax'),
Ns = IntSlider(min=5,max=25 ,step=5,value=5,description='Ns'),
lam = fixed(0.025),
a = fixed(1),
d = fixed(0.25));
interactive(children=(FloatSlider(value=0.0, description='$\\frac{X_p}{a}$', max=2.0), FloatSlider(value=2.0, …
VEGE¶